[f,J] = create_func();
f_obj = [0;0];


results = struct();

k = 1;
epsilon = 0.00001;
max_iter = 100;

% u = [1;0;0]; % initial input
u = [1;2;1];

results(k).u = u;
results(k).u_l1norm = norm(u,1);
results(k).f = f(u(1),u(2),u(3));
results(k).J = J(u(1),u(2),u(3));
results(k).g = results(k).f - f_obj;
results(k).g_l1norm = norm(results(k).g,1);

adaptive = 'true';
% adaptive = 'false';
beta = 0.3;
c_0 = 1.0;
c_max = 2*(1+beta)/(1-beta);
gamma = 0.9;
h = 0.5;

mmt = zeros(size(u));

while results(k).g_l1norm > epsilon
    if k <= max_iter
        
        if strcmp(adaptive,'true')
            c_0 = 1 + (c_max-1)*gamma^k;
        end

        calc_w

        mmt = beta*mmt + (1-beta)*w;
        u = u - 1.0*mmt;

        k = k+1;
        results(k).u = u;
        results(k).u_l1norm = norm(u,1);
        results(k).f = f(u(1),u(2),u(3));
        results(k).J = J(u(1),u(2),u(3));
        results(k).g = results(k).f - f_obj;
        results(k).g_l1norm = norm(results(k).g,1);

        pause(0.01)
    else
        break
    end
end

if results(k).g_l1norm < 0.01
    close all
    find_solution
end